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• ABSTRACT 


A computer program has been prepared to calculate the elec- 
trical characteristics of silicon solar, cells as a function of cell 
parameters and space environmental factors. The program, in For- 
tran IV, computes short-circuit current and the current and power 
at selected voltages to permit construction of a current -voltage 
curve characteristic of the solar cell. Environmental factors con- 
sidered are illumination, intensity and spectrum, cell temperature, 
and exposure to corpuscular radiation. 
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INTBDDUCTIOK 


This is the final report of a program to generate a useful 
mathematical model for silicon solar cell performance in space. 

The model, begun under previous contract*, has been extended, ad- 
justed to a 'closer fit of available measurements, and incorporated 
in 'software for computer simulations of specific designs and ’en- 
vironments for solar cell performance. 

Previous quarterly reports bave presented detail's of the 
work. In this final report, the analytical effort is summarized. 
Beginning with a review of the solar cell and its applications, 
the presentation is arranged to describe first the characteristics 
of the cell (Section 1) , then details of optical considerations 
(Section II) t and then radiation exposure effects (Section III) . 
These are the main parts to the model and computer code. The 
code itself (delivered separately to JPL) is described in Section IV. 
Recommendations for further development are presented in Section V. 


Vc 

JPL Letter Contract 952246, 


under Prime Contract NAS7-1000 
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I. 


MODEL OF THE SILICON SOLAR CELL 


A. Solar Cell Parameters 

Modern solar* cells for spacecraft applications are typified 
by Figure 1. The cell area may be from 1 x 2 cm to 3 x 3 cm, with 
a cell- thickness of 8-14 mils (dimensions are conventionally given 
in these heterogenous units). A junction Is formed about 0.5 micron 
(5 x I0~ B cm) below the sunward surface; this junction separates the 
thin "surface"' region of n-type silicon from the thicker "base" 
region of p-type silicon. For the sunward side, the electrical con- 
tact is In the form of a bar along one edge and a number of thin grid 
lines extending from It across the surface. The dark side contact 
generally completely masks the back surface of the cell. 
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•Surface region (about 0.5 u) 
—Base region (about 10 mils) 


-x 


Si 


Figure 1. 


Exploded view of' typical solar cell. 


Illumination of this device results' in the sunward surface 
being a fraction of a volt positive with respect to the back sur- 
face. A current is thereby induced in an electrical connection 
between the front bar and the baek contact. This is not equal, 
however, to the current generated by the sun (photovoltaic current), 
for some of the photovoltaic current is returned through the solar 
cell itself, in accordance with its diode property. The division 
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of current leads to the solar cell equation, which states that the 
electrical current flowing from a solar cell equals the difference 
between. the photovoltaic current produced in it, and. the diode 
current r D lost in it. 

The output of a solar cell depends on a large number of 'in- 
dependent' factors. First, the intensity of sunlight and the angle 
of incidence are important. A cell directly facing the sun" receives 
the Maximum possible' sunlight,’ If it deviates from this direction 
by an angle 8, as shown in Figure 1, the radiant' energy striking the ' 
cell is reduced by the cosine of 8. As the angle increases, edge effects, 
especially when a covers lide is used, make the reduction deviate 
slightly from this law. Measurement must be 'relied on for each specific 
geometry' when extreme accuracy is required for solar cells illuminated 
at large angles from, the perpendicular. 

The output of solar cells has frequently been characterized by 
the short-circuit current I gc # the open circuit voltage V , and the 
power and voltage of the cell near its maximum power point. While 
these parameters' do characterize the electrical output, in theoretical 
work more fundamental terms are to be preferred. These are- the photo- 
voltaic current , the diode parameters --I and , and the internal 
series resistance R. These four parameters are chosen since it is a more 
straight-forward exercise, -to determine from- semiconductor, physics how 
-they are affected by the environment and 'history of a solar cell. . 

' Observing that a solar cell acts as a diode operating in opposi- 
tion to a current source," Prince and Wolf proposed the solar cell equa- 
tion (ref, 1). Our .only change "to this equation is to replace "the ' 
usual ‘expression AkX/q'by the single parameter V . We propose calling 
V q the ’’characteristic voltage" of the solar cell, and’ write the equa- 
tion for the current I as 


0H-I8) /V 


1 = 1 . 


r \ v-r-jLja./ / v -1 

- I -e °-l. 

O L is 


where V is the voltage across the cell. 


( 1 ) 
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Does the theoretically - inspired Eq. -1 indeed reproduce an 
actual solar cell output? 'We have developed curve - fitting tech- 
niques that select values of the parameters very well. Two samples 
are presented in Figure 2 with the fitted 3olar cell equations-. These 
comparisons demonstrate the accuracy often possible with the solar 
cell equation. They. also demonstrate magnitudes of the parameters. . 
(Voltages here are given in millivolts and currents in milliampereg. 
The resistance of typical cells is generally lower, about 0.1 ohm, 
than in these examples.) Deviations of measurements from the solar, 
cell equation have. been reported. Wolf and Rauschenbach (ref. 3) 
have recommended that the current source be assumed to be shunted 
by two diodes with different characteristic voltages. The revised 
solar cell equation r in our notation* would then take the form 


1 * X L - 


$ 


(v +IE3l )/v ol 




o2 


< ?+ 1E 8 2>/ T 0 2 


-I 


( 2 ) 


This has the immediate advantage of allowing two extra 
parameters for curve fitting, and the disadvantage of relating these 
parameters to the individual solar cell and its environment. The 
rationale for the extra diode term must be developed if it is to be “ 
adopted in a general mathematical model. Recognising that a solar 
cell can have voltage and material gradients along its junction sur- 
face can lead to such composite diode terms. Further, Lad arty; (ref. 4) 
-has derived a silicon diode expression to replace the pure exponential 
in Eq. 1 by a sum of exponentials. These studies, valuable in extend- 
ing our Insight into the action of the solar cell, have not yet led to 
reliable formulas for engineering use. Although a higher accuracy in 
matching a given I-V curve can be achieved, the simpler solar-cell 
equation with its four parameters is to be preferred for parametric 
studies . 

The space environment operates on the solar cell parameters in 
ways that are imperfectly understood. Of course, the uncertainties 


4 




Figure 2. 


I « 60.2 - 
I =* 46.3 - 


Comparison of measured I-V .curves (rpfi .2). with 
points computed with the solar cell equations : 


1.675 x 10-*(e ( v+ *- 791I)/45.34 -]L) • ceir c5c _ 66 

1.625 x 10~ 4 (e (▼+ -3220/46.98^ Cftll T5 _ 6 
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in radiation fluences in space add to the difficulty. As a result, 
theory and experimental data from the laboratory are often the best 
that is available to predict changes in solar cell output in space. 

By correlating theory with laboratory data, the relation of these 
four parameters to environment and exposure can be obtained. 

B. Photovoltaic Current Density 

To determine the magnitude of the photovoltaic current produced 
by the solar cell, it is necessary to follow the sequence of events 
from light absorption, through carrier diffusion and drift, to recom- 
bination at the junction. Damage due to corpuscular radiation and the 
dependence of the physical parameters on temperature and illumination 
are to be taken into account in this sequence. 

The equations obeyed are the continuity equation for the minority 
carrier current 

i m . Hisx + G(x) .0 <3) 

q dx . ?(x) 
the current equation 

J =* q D~| + q p,E n (4) 


the damage equation 


L? 


+ K § 


(5) 


* 


the symbols used throughout are given in the glossary; 
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the diffusion relation 

L 2 = TD (6) 

« H", 

and the Einstein relation. 

q B = |i kT (7) 

These equations, applicable to the minority carriers in the base 
region, have their counterparts for the surface region, using a sepa- 
rate set of values for the material parameters. The photovoltaic cur- 
rent I is, of course, the sum of the current densities J into the 
junction from both sides., multiplied by the area of the solar cell that 

is exposed to the light source. This is somewhat less than the actual 
■ ‘ * < 
front surface area, due to masking of up to 107c. of the surface by the 

front' contact. 

The effect of temperature on the diffusion coefficient D of minority 
carriers in silicon may be determined 1 from measured values of the temper- 
ature-mobility relationship. The measurements (ref. 5) show that for 
electrons the mobility varies as T~ E * 5 , and "that for holes it varies ag 

T" 3 * 7 . Therefore, D varies as T 1 ‘‘ 5 and D varies as T“ 1>7 . 

- , n . - p 

The effect ,o,f temperature of the minority., carrier lifetime t is 
composed of .many factors, most of which depend- upon the. specific pro- 
cessing that the cell has been subjected to.. Each of the process steps 

«... . • " * ' ■ x *’ 

changes .the number,, type., and energy .level of the recombination centers 

i * r 

that determine the lifetime. The temperature dependence of many types 
of recombination centers has not been described in the literature. It 
was therefore necessary to assume a relationship that is consistent 
with experimental results. If -r is taken as proportional to T S,S , the 
temperature coefficient of I in the region about 300°K is 0.057 o /°C, 
which is consistent with the photovoltaic current measurements given 
by Reynard (ref. 6) and others. This back - calculation of temper- 
ature - dependence of t from observed temperature - dependence of 1^ is 
an empiricism worthy of further review. 
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1. Base Region Contribution 


The current equation is simplified in the base region since 
the electric field E is negligible in a uniformly doped region.' 
Combining the current equation and the continuity equation then 
yields 


D “ + G < x >' * 0 < 8 > 

n dx 2 T 

We assume for boundary conditions that the minority carrier density 
n vanishes at the junction and at the cell surface. Equation 8 is 
then solved by computer iteration using a variable interval for the 
distance x in a difference approximation to the differential. 

A variable Increment technique is used to increase the accuracy 
of the difference approximation and to minimize computer time. The 

* '• ‘5 * 

difference equation approaches the differential equation as the in- 
terval size approaches zero. Therefore small increments, particularly 
near the junction where -the current equation (Eq. 4) is used, are de- 
sirable. However, use of very small increments throughout the solar 
cell thickness requires a large number of increments and, consequently, 
more computer time. For this reason a compromise was made: small in- 
tervals in the critical region near the junction and increasingly larger 
ones in the less significant areas deeper, into the cell. The technique 
necessarily requires a more complex form of the difference equation 
-than normally used. 
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n 



Figure 3. Construction of. variable mesh to approximate 
the curve n(x) at points 


The difference approximation, for a variable width interval, 
can be derived with the aid of Figure 3. The slopes of the curve 
between a' ’'and b, and between b and c are approximated as constants. 

The first derivative of the curve, evaluated at b, is taken as the 
average of the slope on either side of b, i.e. 

*5 I I ”kH '"k N i _ °!Hr2 h k^ n fcflV bri-l ' bt )~°kV:l ' 

\h \ > 2h kVl (9) 

l * ’ ' 

The second ..derivative at b is determined from the slopes on 
either side -of b. ^..It is then 

A 3 n [ „ \ " n k-H * h k-blJ +n k \n nr 

ts c 3 b \ h k+l [\ * '\+ly " T V 

Using equation 9 and 10 in the continuity equation, with the 
field terms eliminated, the minority carrier concentration can be 
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found from 



To solve Eq. 11, it is necessary to have values for and 
nj . We assume the boundary conditions that the carrier density 
vanishes at the junction and at the cell surface. This assumption 
was also made by Bullis and Runyan (ref. 7) but other boundary 
conditions have been assumed. We guess n x and calculate all the 


higher values of n^. 


We repeat this guessing of n^ until we 


arrive at a satisfactory value for the carrier density at the cell 


surface. 


The accuracy of the initial n^ and the iteration technique 
are of importance in determining how often the calculation must be 
repeated before obtaining zero carrier density at the contact sur- 
face of the cell. If the guess fornj. is too small, then the values 
of n fc determined via Eq. 11 will change sign in the cell. If it is 
too large, the at the back of the cell will fail to be zero. 
Iterating on n^ leads to as close an estimate as is desired. One ■■ 
possible technique for convergence is to compare each with ^ 
and if there is a sign change then stop, increase the estimate for 
Uj by a nominal 10% and repeat. When there is no sign change, de- 

a 

crease the estimate for n, by a nominal 5% and repeat until a sign 
change occurs. Then increase by 1% until there is no sign change. 
Such a convergence routine can obviously be carried to any level cf- 
accuracy in the estimate of vt^ for the solar cell In question, by 
taking advantage of this sign change. Our present technique is 
described on page 31. 
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The current, equation put into this difference form relates 
nj »• the minority, carrier concentration at the distance from 
the junction, to the current from the base into the junction. 

. J = q D n /h, (12) 
H n x l 


since the boundary conditions require that n vanish at the junction. 


• 2.*, . Surface Region .Contribution ; 

The surface region of a solar cell, typically only 0.-2 to 
0.5 microns' thick, contributes a minor amount to the photovoltaic 
current.' Calculation of this contribution with the technique- 
used for -the base region is- complicated by 'several factors. To 
begin with, "the normal process- of junction formation by in-diffus- 
ion of n-type impurity ‘atoms results in a complementary error 

, » 

function type- of' distribution of donor atoms. The concentration 
N(x) is related to the concentration N(0) at the surface by the 
t emp era ture -dependent ''Et&jflic diffusion coefficient p and the pro- 
cess time t by 

H(x)'=N(0) erfc x/V^t J - (13) 

» . s 

In solar cell manufacture, a typical value of N.(0) appears to .-be 
10 2 ° atoms/cm 3 . The product pt, if unknown from the- process, can 
be calculated if the junction depth is known. At the junction, 

N(x) must equal the impurity concentration' of the base material. 

The' code uses a polynomial approximation to the complementary 
error function. 


erfc 


■ 60 - ( 


c nf + chi 3 + c n a - + c. n) 
4 3.3 l / 


2 

VTT 


- y 


(14) 
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where 


. 1 

n = 1 + 0.381965 y 


c t = 0.12771538 
Cg — 0.54107939 
c 3 = 0.53859539 
c 4 = 0.75602755 
y ~ x/V4jbt 


The non-uniform impurity concentration over the surface region 

also causes variation in the minority carrier diffusion coefficient 

D . The term D is the diffusion coefficient for holes and should 
P P 

not be confused the atomic diffusion coefficient used above. 

Since does not change uniformly with concentration, the variation 
with impurity concentration shown by Gonsrell (ref. 8) is .used. 

Since no simple equation will fit the entire curve, it was decided 
that the best approach was a series of empirical equations. These 
are: 


N 

< 10 1S 

D 

«s 

13.0 

(a) 



P 



10 1B 

< N < 10 16 

B 

P 

=s 

52.0-2.61og N 

(b) 

10 16 

< N < 10 17 

d' 

=3 

56.16-2.86 log K 

(c) (15) 



P 


- 

10 17 

< N < 10 18 

D 

P 

sa 

82.86-4.42 log N 

(d) 

10 18 

< N < 10 19 

D 

P 

tx 

41.28-2.1 log N 

(e)" 


N < 10 19 

D 

T* 

— 

1.0 

(f) 


For the surface region, where the field is not negligible, 
the continuity equation becomes 


G(x) 


i. ■ i 

T P 


P H- 


dE 


p dx 


- u E -y 4 - + 
■ 'p- dx 


P dx 3 


(16) 
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where P = concentration of minority carriers (holes) 

T ** lifetime of holes 
P 

mobility of holes 

• Vm .1 Sh 

E = field due to impurity gradient = ~ — ~ 

Using equations 13, 14, 15 and 16 we solve the continuity 
equation on the computer by putting it in- the form of a difference 
equation, in the manner demonstrated in the previous section. 

The inclusion of- these equations in the final program would 
unnecessarily lengthen the computing time.. The routines for the 
surface region, therefore, were run separately, with varying 
temperature, junction depth, radiation and illumination intensity. 
The variation of the surface contribution to the minority current 
density J with these factors was curve fitted as 

J = (0.01408 + 0.005956 0n - 0.01411e “ 5 - 383 ^ 

(T + 400) (U cos 9/ 140) mA/ cm 2 

The minority carrier diffusion length L (in microns) is taken at 
the first 'increment in the base region. The assumption here is 
that L, computed in the base region near the junction for base 
region calculations, is a measure of the amount of radiation ex- 
posure the surface region has received. For the thickness of 
practical solar cell surface regions, the damage is uniform and 
this value of L is an indicator of the exposure. 

Use of this expression to estimate the surface region con- 
tribution to the photovoltaic current permits us to bypass a 
more involved calculation such as used for the major contribution 
by the base region. 
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G. • Diode Characteristics 


The characteristic voltage of the diode response by the 

solar cell appears to be- independent of temperature, radiation 

exposure, and illumination. For this reason', the expression AkT/q 

normally found in the solar cell equation has been replaced by a' 

constant V in our solar cell equation (Eq. 1) . This implies that 
o - 

A is inversely proportional to temperature, and measurements by 
Kennerud (ref. 9) replotted in Figure 4, back this up. The lack 
of any clear dependence of V q on radiation exposure was demonstrated 
earlier (ref. 10) . Analysis of solar cell response curves under 
different illumination intensities has also failed to show any clear 
changes in V . 



T( e K) 

Figure 4. Reported values of A, versus temp- 
erature T , for typical n/p silicon 
solar cells. After reference 9. 
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The antithetic situation occurs with respect to the diode 
saturation current I . A functional relationship to light intensity^ 
temperature » and radiation exposure has been developed in this study 
from a combination of theory and indirect measurements. 

Relating the diode saturation current to illumination can be 
accomplished from experimental data on the open circuit- voltage -V - c> 
remembering that the photovoltaic current is proportional to illumi- 
nation. According to Ritchie and Sandstrom (ref. 11) V increases 
at a rate of 0.2 mV/mW/cm 3 with increasing illumination intensity IJ. 

I can be calculated from this using the standard solar cell equation, 
o * - 

If I is plotted against 9m 0, as in Figure 5, the curve can be approx! 
mated by a straight line, tie therefore make 1^ proportional to 9nV 
in the mathematical model. 



Figure 5. Variation’ of diode current 

with illumination intensity. 
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A similar technique was employed to establish the temperature 
variation of I . Theoretically I is proportional to exp (-E /2kT) . 
Assuming the photovoltaic current density increases $,057o/ G, and 
the series resistance is negligible , the theoretical relationship 
yields a temperature coefficient of V as 0.56%/°C. This value 
falls within the range of Reynard’s early experiments (ref. 6) in- 
dicating the correctness of the theoretical expression. ’ ' 

Shockley’s analysis (ref, 12) for diode junctions provides a 
theoretical expression for I that is inversely proportional to the 
base region minority carrier diffusion length near the junction. 

This establishes a radiation dependence of 1^ through the damage 
equation (Eq. 5). Thus, we expect the diode current of a solar cell 
to increase with increasing radiation damage. This trend has been 
observed (ref 13) in an analysis of ATS-1 spacecraft solar cells. 
Coupling this with the temperature and illumination intensity relat- 
ionship yields, in mA/cra 3 

X o * <L 0 /L)[1.57 fci (U cos 6) ~ 3.36] exp (-6492/T)x 10 B (18) 

The factor cos 9 relates the solar illumination intensity U to the 
intensity observed on the solar cell surface as a function of the 
aspect angle 9. This expression assumes the cell is illuminated on 
its front surface. hare is a current density and is to be scaled 
appropriately to the particular solar cell area under consideration. 

D. Cell Resistance 

The series resistance of a solar cell is a composite of resis- 
tance terras due to current flow across the bulk region, along the 
surface region t© the contacts, and from silicon to front and back 
contacts (ref. 14). Each of these terms logically would behave 
differently under radiation exposure. As a consequence, careful 
analysis of resistance effects would require one -first to partition 
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R among its components, and then scale these components as they 
change with radiation exposure. Some uncertainty is obvious in 
this scheme (even if the effect of radiation on silicon resistiv- 
ity were known) , since manufacturing tolerances would vary the 
components from cell to cell. Finally, discontinuous events such, 
as the lifting of a contact would be difficult to predict. This 
would lead to a value of R much greater than normally expected. 

Cell resistance generally ' plays a small part in the behavior 
of a solar cell. In this study it is considered a constant due to 
lack of information on its variation with radiation exposure. Ex- 
perimental measurements would be useful in establishing the quan- 
titative nature of the increase in R. Whether the nature of the 
dopant affects the change, whether it is linear with exposure, 
whether temperature plays a significant role, and whether the 
nature .of-, the damage depends on the bombarding particles, are quest 
ions to be answered by such measurements. 
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XI. CARRIER GENERATION 


The light generated rate of production of minority carriers 
G(x), occurring in the continuity equation (Eq. 3), is a function 
dependent on the illumination, intensity', spectrum, and absorption 
coefficient of light in silicon. The non-analytic variation of 
spectrum and absorption coefficients with wavelength of light pre- 
cludes an analytic solution. Therefore, the contributions to G(x) 
from different parts of the light spectrum must be considered in- 
dividually. This is accomplished through numerical Integration of 
the relationship 

g<*. x> - «(1 )h«|h!) e - ® * /c ° se ' 


over all values of wavelength X to which the solar cell responds. 
This ‘ region is normally taken as 0.4 to 1.1 microns. The small 
fraction of incident light reflected at the solar cell surface, 
normally on the order of two or three percent, has been neglected 
in this expression. 

Values of the absorption coefficient or(X) and the spectral 
irradiance H(X) are given in Table 1 for an illumination intensity 
of 140 mN / cm 3 . AMO (air mass zero) conditions, and a temperature 
of 300°K, Variation of the illumination intensity causes a pro- 
portional variation in H(X) for all X. Variation in the air mass 
condition, e.g., reduction of the spectrum to AMI conditions by 
atmospheric absorption, varies the spectrum since the absorption is 
selective and changes with atmospheric conditions. In this case the 
individual values of H(X) must he changed appropriately. Temper^ 
ature will effect the absorption coefficient. Macfarlane and 
Roberts (ref. 17) give the absorption coefficient as 


1 

f fuu - E - kB 
£ 


f litu - E -f kB ") 
g 

2- 

[l-e- B/T 1 

^ lieu 

j 

j 
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of = A 


( 20 ) 



Table 1 


Absorption Coefficient of Silicon and Sunlight Intensity 
as a Function of Wavelength 

k A, in microns, a ( X )' in cm -1 (Ref. 15) 

R(Jl) in watts/ era 2 - p, (Ref. 16) 


X 

■Q'(X) 

H(X) 

0.40 

7.50 x 10 4 

0.1540 

0.45 

2.58 x 10 4 

0.2200 

0.50 

1.18 x 10 4 

0.1980 

0.55 . 

7.00 x 10 s 

0.1950 

0.60 

4.65 x 10 3 

. 0.1810 

0.65 

3.33 x 10 s - 

0.1620 

0.70 

2.42 x 10 s 

0.1440 

0.75 

1.69 x 10 s 

0.1270 

0.80 

1.12 x 10 s 

0.1127 

0.85 

7.95 x 10 s 

0.1003 

0.90 

3.80 x 10 2 

0.0895 

0.95 

1.80..X 10 2 

0.0803 

1.00 

7.30 x 10 1 

0.0-725 

1.05 

2.08 x 10 1 

0.0665 

1.10 

4.40 x 10° 

0.0606 
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where : 


A is a constant 
B - 600°K (for silicon) 

"fito — energy of incident light (in eV) 

E = energy gap of the semiconductor (1.1 eV. in silicon) 
S 

k = Boltzmann s s constant 


Since E^ varies very slowly with temperature (-0.00013 eV/ C for 

— j 2 

(fits) - E - kS) /fits | and 
-.2 6 J 


silicon) (ref. 18) the quantities 


[ 


(ftcu-E 4- kB) /tuo 


can be considered constant' for each wavelength, 

At X - 1.1 \i , 


making or proportional to ; l/(l-e ~ + l/(e B ^ T -l) j 


where fits is very close to E^, the value of a changes only 7% between 
250°K and 300°K. Since light in this region of the spectrum con- 
tributes only a small number of carriers, the effect may be neglected. 
At lower wavelengths the variation is even smaller. If the above quan- 
tity is designated as f , and (f - 1) is plotted as a function of B/T, 
as in Figure 6, a straight line results. The equation which fits this 
line with a maximum deviation of leas than 1%, is 


bn (f - 1) * 1.201 - 1.171 


J 1.6 &ZA 

T i 


( 21 ) 


From' equation 21, the relationship between ca and T can be shown to be 


i -702. 6/T,- 

a (X,T) = a (\,300) — (-22) 



Figure 6. 


Relationship of absorption 
coefficient to temperature. 
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III. RADIATION EFFECTS 


A. Proton Shielding 

The thickness of a solar cell is comparable to the distance a 
proton can travel in silicon when its energy is typical o£ protons 
found in space. Shielding by a coverslide and self-shielding 'bjr the 
solar cell are consequently important. This is especially signifi- 
cant when one considers the proton energy - dependence of damage. 

s 

The distance of travel, or range R, is often related to the in- 
cident proton energy E by formulas of the form 

“ : * 

R « R E n (23> 

o o 

The equatioh is not exact, but good fits can be provided over limited 
ranges of E . -Table 2 is such a fit to tabulated data. 


Table 2 

Values of R q and n for proton range - energy relationships and corre- 

j > • •& 

spending energy intervals. 


Energy (MeV) 

R o (mg/cm a ) # 

n 

0.0 £ E < 0.3 

2.81 

0.995 

0.3 £ E < 0.8 

3.95 

1.277 

0.8 ^ e < 2.0 

4.11 

1.460 

'-2.0 5 E -< 200 

3.42 

1.726 


computed from data in ref. 19' 
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The range formula implies that the energy of a proton along its 
track can be calculated from the residual distance it is to travel 
before stopping. The relation is not exact, for there is some strag- 
gling of the individual tracks of protons of the same energy, but this 
is generally quite small. The average .straggling, as a fraction of R, 
decreases with proton energy E q and it is lass than 4 % for 100 keV 
protons (ref. 19). Thus, we treat the range equation as being exact 
and compute a proton energy (E n - — } / for protons of initial energy 

E which have traveled a distance x in silicon. When monoenergetic 
o 

protons of an omnidirectional fluence strike the solar cell surface, 
the effects of slant penetration cause a spectrum at depth x. 

A covers lide of thickness t will remove protons of energy E and 

Q 

incident angle 8 with the normal if their range R is less than the 
path length t/cos 8 through the coverslide, .It will also reduce the 
energy of a transmitted proton to E given by 

E » j*E n “ RcosT |^ n (24) 

L. O -i 


The differential proton spectrum striking the solar cell due to 
a monoenergetic, isotropic unit flux penetrating the covers lide can 
be derived from thi3 expression and is 


§ 

P 


(E) = 


ntE 


n - 1 




R. 


(25) 


for energies E less than the energy of a proton penetrating at normal 
incidence, as calculated fromEq.. 24. 

The total proton fluence at any depth z in the solar cell 
assembly is related to the proton fluence in space by 


§ 

P 


(E o ,2) 


* nz 


A e (E )E n_1 

i — 

Jo R (E n -E*) 2 
o s o 


d£ 


. (26) 
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These equations permit calculation of localized proton damage 
in solar cells with covers lides, given the appropriate damage 
equivalence. This is presented in the next section. 

B. Proton Damage 

The proton damage coefficient is the measure of decrease in 

minority carrier diffusion length due to a fluence § of protons (K^ 
equals the incremental increase in the quantity 1/L 3 with incre- 
mental Increase in §) . The damage is due to Rutherford scattering 
of protons, which dislodges silicon atoms from their lattice position. 
Hence, the proton energy dependence of is approximately given by 
1/E, which is the energy dependence of the Rutherford scattering cross 
section. The threshold for dislodging atoms, corresponds to a minimum 
proton energy of about 0.0001 MeV. 

Crowther, et al. (ref. 15) have found a flattening of the energy 
dependence of K below about 0.5 MeV. This effect, may correspond to 
an annealing mechanism whereby a dislodged atom has not been pushed 
far from its site, and has a high probability of return. The follow- 
ing equations fit these measurements for 1 ohis-cm p-silicon and pro- 
vide a ratio for higher resistivity p-si-licon that agrees with measure- 
ments by Denney and Downing (ref. 20) . 

K (E) = 1.2Q~°-' 7B E “ °* 9 x 10“ s , [e > 3 MeV 1 - (27) 

P . , . *- . J 

K^(E) =1.92 Q~ 0, ' 7B e * 1 ' 04 (E/.962) ~°* sB s: 10“ S [3>E>l] (28) 

K p (E)=U92TI“° -7B e “ 1 * oaE x 10“ B jjL>E>10~ 4 j (29) 

Because proton energy changes rapidly with depth of penetration 
into the cell, and because the damage coefficient is so dependent 
on proton energy, damage by protons of less than about 5 MeV results 
in a ' highiy notiuniform minority carrier diffusion length across the cell. 


23 . 



C« Electron Shielding 

When a solar cell is covered with t grame/cm 3 of coverglass 
the electron f lux at the cell surface due to isotropically incident 
moaoenergetie electrons of energy E is a spectrumover lesser energies 
and is dependent on the thickness of the shielding. Here* however, 
deflections in the individual electron paths complicate the analytic 
determination of the shielding effect, k weighting factor, derived 
from a study of Mpnfce Carlo results (ref. 21)* enables us to deter- 
mine the effective damage coefficient IC e (E,x) at depth x into the 
cell as a fraction of the damage coefficient of the incident elec- 
trons K e (E) . The weighting relationship is given by 

K (E,x) = K (E) exp [ - I0(x + t) /E 1 * B ] (30) 

C £ 

where xl t Is the areal density in gm/cst; of material penetrated to 
the site of the damage being considered.' 

D, Electron Damage 

The damage coefficient K £ f©r electron damage in p-type silicon 
has been fitted empirically as shewn in Figure 7. We neglect to 
study fits for n-type silicon since electron damage effects are neg- 
ligible in the surface region, compared to those in the base region. 

We plotted^h^ square root of measured values of K versus electron 

£ 

energy. For crucible - grown p-type silicon of resistivity 10.6 ohm- cm, 
the measured points, shown in Figure 7, can be connected by two 
straight line segments. That a straight line results over the energy 
1-40 MeV indicates that a recombination center requiring two defects 
may be involved (ref. 21). Below 1 MeV, sufficient data are not avail- 
able for such a conclusion, but a straight line curve fit can be pre- 
sented. 

This fitting of the data from ref. 22 and the dependence on elec- 
trical resistivity discussed above, give the damage coefficient 8% 
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K (E) = (10/fl)°* 6 (1.2 Pm 2 . 17E) s x 10 ' 

E > 1 MeV "| 

(31) 

K (E> = (10/0) 0,6 (0.67&J 4E)' s x 10~ K ! 

u ! 

1 >E > *25 j 

(32) 

K e (B) - 0 | 

i .25 >E j 

(33) 


where 0 is the resistivity in ohm-centimeters and E is the electron 
energy in MaV. The expression to fit measurements between 1 and 40 
MeV will underestimate the damage for lower energies. Below 1 MeV, 
the second factor in equation 31 can be replaced by (0.67 9n 4E) 2 . 
This fits the measurement at 0.6 MeV and the generally observed 
’'apparent" threshold of 250 kaV. 
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IV. PROGRAM MODULE DESCRIPTIONS ' 

The theory described in the previous sections has been in- 
corporated in a computer program to provide estimates of the 
simultaneous effects of environmental factors on solar cells. 

The function of the program la twofold: the. determination of the 
solar cell parameters and and the . associated I-V curve 
under, various environmental conditions, including radiation 
damage by electron and proton spectra. The program, is modular, 
i.e, , a main program accompanied by several subroutines. The 
main program collects the data necessary to describe the solar 
cell and its > environment, and calls each appropriate subroutine 
when needed. Each subroutine performs a specific- task and is 
discussed individually below. 

The environmental input consists of the temperature, the 
illumination intensity, and up to ten each of proton and electron 
energies with their associated fluenees. Each fluence is con- 
sidered to be incident isotropically. , By suitable choices of 
energies and fluenees, continuous particle spectra can be approxi- 
mated.' Particle energies are restricted to values below ZOO MeV 
for protons* and below 40 MeV for electrons. Particles with energies 
above these limits are ignored by the program due to lack, of data 
on high energy damage coefficients, ’ r . ' 

Optical parameters, as given in Table 1, are stored in the 
program. A list of other necessary Input variables and- their 
assumed values is given as Table 3. The values may be changed as 
desired, but care must be taken to retain the specified dimensions. 
More detailed discussions follow in the subroutine descriptions. 
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Table 3. Input Variables 


Variable 

XLO 

D 

RHQ 

VO 

R 

T 

XJU 

CT 

TEMP 

U 

THETA 

EOP (I) 
EHP (I) 
EOE (I) 
PHE (I) 

M 

INA 

A(I) 

H(I) 

GIN 


*(I in this 


Definition 


initial base region minority carrier 
diffusion length at 300°K 

base region minority carrier diffus- 
ion. coefficient at 30Q°K 

base region resistivity 

.solar cell characteristic voltage 

solar cell series resistance 

solar cell thickness 

junction depth ‘ ’ 

covers lide thickness 

temperature 

illumination intensity 

incident angle of illumination with 
respect to normal incidence 

proton energies 

proton fluences 

electron energies 

electron fluences 

- ■_ ‘IX ,• j"' ' 

number of points considered in the 
difference equation solution 

number of beams approximating- each 
isotropic proton fluence t- 

absorption coefficients of light in 
silicon 

spectral irradlance obtained from 
the Johnson spectrum 

index of refraction .of light in 
silicon 


table is a running index,) 


S-tored Value 

{ 150 microns 

35 cm 3 /sec 

10 ohm-cm 
43 rav 
0 . 1 ohm 
14 mils 
0.5 microns 

i 

6 mils 
300°K 
140 mW/ cm 2 
<* 0 . 0 ° 

0.0 MeV , 

0.0 protons/cm 3 
0.0 MeV 

0.0 electrons/cm 3 
200 

» i'i, . > 

50 ' V 

* * J « “ * 

see Table 1 

-■ 'i f *- 

see Table 1 
3.8 ■ 
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A. Subroutine ABSCIS 


Subroutine ABSCIS computes the increment thickness HX(K) 
necessary for the difference equation calculation, and the depth into 
the cell EX(K) of each corresponding point. The first 20 increments 
are of equal width delta; the remaining values of HX(K) are given by 
(K-20) delta. Therefore any error induced by the unequal increment 
technique occurs away from the junction. Such a compromise between 
equal and unequal increments reduces the degree of approximation in 

f'5 ~ ,w 

the critical region near the junction, resulting in a more accurate 

) * 

evaluation of the minority carrier concentration and, consequently, 
the photovoltaic current density. 

The number of points M is set at 200. Testing the program with 
a greater number of points (smaller increments) shows an insignifi- 
cant change in the output. 

B. Subroutine LIGHT 

• SJ 

The rate of production of minority carriers per cm per 
second due to light absorption is computed for each point provided 
by subroutine ABSCIS . A Simpson’s rule, integration over the AMD . 
Johnson spectrum (ref. 16) from 0.4 to 1.1 microns in 0.05 micron" 
steps i, s. employed. Other spectra, e.g., the tungsten spectrum, may 
be used by. changing the spectral lrradianees H.(I) in the subroutine’s 
input data. to those of the desired spectrum for the wavelengths 0.4, 

0.45, 1.1 -microns. The method assumes that in any two 

successive 0.Q5. micron steps the spectrum approximates a quadratic 
function. Where this is not the case, e.g., the xenon spectrum, 
smaller intervals must he employed, or the spectral pead^jl- averaged 
in the Interval. 

C. Subroutine COYER 

COVES, approximates each monoenergetic isotropic proton 
f luence by a a.Ct of ,50 beams incident at angles ranging from.. zero. 
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with respect to normal incidence, to the maximum angle a proton of 
the given energy can have and still penetrate the coverslide and 
surface region. This yields a maximum angular' increment of 1.8°. 

Finer increments can be used but the small change in the output does 
not justify the additional computer time. For each angle the proton 
energy after penetration and the incremental fluence over the asso- 
ciated angular increment are determined. Computation terminates if 
the initial proton energy is insufficient to penetrate at normal 
incidence, or if the energy after penetration is insufficient to cause 
damage. 


D« Subroutines PROTON and DAMAGE 

Associated with each of the 50 beams determined by COVER is 
an angle dependent proton energy profile through the solar cell thick- 
ness. PROTON utilises the range - energy calculations of Janni (ref. 19) 
to determine the energy at each point in the cell until such time as: 1) 
insufficient energy remains at a point HX(K) to penetrate the next in- 
crement HX(K) ; 2) the energy is below the damage threshold; or 3) the 
cell has been completely penetrated. ‘ It is assumed that the proton 
follows a straight path until one of the above conditions is met. To 
each of these energies subroutine DAMAGE associates a damage coef- 
ficient, These, in turn, are used with the incremental fluence to de- 
termine the degraded minority carrier diffusion length. A damage 
diffusion length profile as a function of position in the cell results. 
Since the process is repeated for each beam of each initial energy the 
final diffusion length profile represents the total damage done by the 
approximated proton spectrum. 

E. Subrouting EIiECT 


ELECT is the electron counterpart to subroutines COVER* PROTON , 
and DAMAGE.' . Here, however, the assumption of a non-doflected path can- 
not be made. We instead use a weighted damage- coefficient to represent 
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the damage at a given depth EX(K) in the cell* as described on 
page 24. The damage coefficient of. the incident monoenergetic 
electron -finance is- computed first., .The effective damage coef- 
ficient at depth M(K) is then determined as a function of the 
areal density in gm/cm 2 of material penetrated, i.e., EK(K) plus 
the coversll'de thickness. The coverslide and silicon densities 
are 2.2 and 2.33 gra/e!a% respectively. 

As in subroutine DAMAGE the minority carrier diffusion length 
profile is updated using the associated electron fluence and the 
weighted damage coefficients. The process is repeated for each 
electron energy given as input. 

F. Subroutine BOOT 

Subroutine SOOT contains the iteration technique to solve 
the difference equation approximation of the continuity equation. 

The value of the minority carrier concentration at the junction, C(l) * 
is set equal to aero and an approximation of 10 s made for its value 
at the second point c(2) . The approximation of C(2) ia then increased 
or decreased by an order of magnitude depending on the sign of the 
carrier concentration as computed for the back surface of the ceil 
with the difference equation. A positive sign results when C(2) is 
too .large and a negative sign when it is too -small. . The process con- 
tinues until the signs differ for two success ive r approximations* in- 
dicating that the correct value of C(2) lies between them. The bi- 
section technique is then employed until either the carrier concen- 
tration at the back of the cell becomes zero* or two successive values 
of C(2) differ by less than 10” 4 percent. 

G. Subroutine TRAP 

TRAP contains the actual difference equation approximation of 
■the continuity equation. It.iSicalMSbJ' subroutine KGOT and control 
is returned if: 1) for a particular C(2), the minority carrier 
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concentration anywhere in the cell exceeds I0 ao carriers /cm 3 ; 2) 
the minority carrier concentration becomes negative? or 3) the 
computation for every incremental point is completed. 

H. Subroutine CURVE 

CURVE utilizes the previously determined data along with 
the input data to generate the photovoltaic current density I , 
the diode saturation current 1^, points on the resultant I-V curve, 
and their associated power P. The base region contribution to the 
photovoltaic current density is determined from the minority carrier 

concentration gradient at the junction, C(2) /HX(1) . The surface 

\ 

region contribution is calculated from the junction depth, illumination 
intensity, temperature, and degraded minority carrier diffusion length 
at the junction. Addition of the two current components yields the 
total photovoltaic current density. Temperature, illumination in- 
tensity, and degraded minority carrier diffusion length at the junction 
are used to approximate the diode saturation current. 

A set of the four solar cell parameters is now available. The 
program is completed with the output of these parameters, environ- 
mental data, points on the resultant I-V curve computed from the solar 
cell equation j and the power associated with each point. An example 
of the output is given in Table 4. 
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Table 4. Program Output (Sample) 


I-L = 

40.45 MA 

U= 140 

.0 m/cm* 2 

VO - 

43.00 MV 

TEMP= 

300'. DEG* K 

R = 

0.100 OHM 

THETA 0 

0.00 BEG. 

10 « 

0.76E-03 MA 




I(MA) 

V-(MV) 


0.00 

530.9 

.0.00 

4.04 

■525.9 

2.13 

8.09 

520.5 

4.21 

12.13 

514.3 

6.24 

16.18 

507.3 

s; 21 

20,22 

499.1 

10.09- 

24.27 

489.1 

11.87 

28.31 

476.3 

i3"49 

32.36 

458,4 

14:83 

36.40 

428,2 

15.59 

36.81 

423.7 

15.59 

37,21 

418.6 

15.58 

37.62 

412.8 

15.53 

38.02 

.406.1 

15.44 

38.43 

398.2 

15.30 

38.83 

388.6 

15. '09 

39.23 

376.2 

14;76 

39.64 

358.7 

14.22 

40,04 

328.9 

13-.17 

40.45 

1.1 

0.05 


The information presented here are the four parameters of the 
solar cell equations given a light - intensity U of 140 milliwatts/ cm 3 , 
a sun orientation THETA of 0.00 degrees from the normal* .and a 
temperature TEMP of 300 K. TSfe'jparamefcers are used by the program to 
compute voltage V and power P for twenty values of the current X 
through the external load. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


A. Conclusions 

Wo cow have detailed mathematics available to study solar cell 
capabilities and design. It transcends previously utilized techniques 
of Vowing the solar cell In terms of equivalences * such as represent- 
ing non-uniform proton damage as equivalent 1 MeV electron damage, 
and allows analysis on the level and in the scope of the particular 
phenomena of interest. , For example, the model is capable of predicting 
the power degradation of solar cells in a combined environment of 
protons and electrons both of arbitrary energies. Nonuniformity of 
damage is permitted in the calculations. In the process, it provides 
the necessary framework to optimize a solar cell coverglass assembly 

particular radiation environment, taking into account the other 
environmental factors' and cell parameters. One measure of the versa- 
tility of the model is the list of variables. Table 3, considered by 
the program. 

The computer program encompassing the mathematical model has 
been used to test the validity of the model relationships. The results 
are, in general, consistent with experimental evidence as illustrated 
In previous reports under this contract. Summarizing the results of 
the testing process wa'find: 

1. V oc decreases with increasing temperature, primarily due to 
the increase in diode saturation current. 

2. P decreases with increasing temperature, illustrating that 

max 

thermal effects on this portion of the solar cell I-V curve 
are dominated by I rather than I T . ' 

O -u 

3. I ,7 , and ? decrease with increasing junction depth 

S C OC 1S33C 

beyond 0.2 microns. 

4 . The surface region is essentially unaffected by irradiation. 

In the course of this contract many improvements have been made 
•.in the mathematical model generated under JPL contract 952246. 
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Noteworthy among these are: 

X. A technique of unequal increments has been employed in 
the difference equation calculations, which reduces the 
degree of approximation in the critical region near- the 
junction. This results in a more accurate evaluation "of 
the minority carrier concentration and, consequently, the 
photovoltaic .current density, 

2. The parameters aascciatedwith the surface layer have been 
included, reducing the degree of approximation- formerly 
introduced by assuming a. constant surface region contri- 
bution to the total photovoltaic current, 

3. The analysis of the electron shielding by a coverslide 

has been greatly simplified, facilitating its incorporation 
into' the computer program without appreciably extending 
computation time, 

4. The study of the diode reverse saturation current has 
lead to a more detailed expression that takes account of 
the response of this ‘parameter to changes In temperature, 
illumination, and radiation exposure. 

B. 'Recommendations 

1, Further refinement of the model produced in this work is 
recommended. Experimental efforts tailored to an , investigation -of 
specific parameters -are needed. When one parameter is being 

• i ? - - 

changed, the others should be measured so that their values can 
be -placed in the computer calculation. (The l 3 ck t of^e«apIete t ^ 
measurements makes it difficult to assess many published experi- 
mental data.) An experimental program- would not only validate, 
the model, but would point out areas where future analytical effort 
is warranted.- •- 

2. With respect to 1 sw energy -protons we- have found that 

the theoretical the actual damage to the 
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photovoltaic current in the cases we have analyzed. As discussed on 
p. 23, the energy dependence of the damage coefficient diminishes at 
low energies (below 1 MeV) . A reasonable extension of this behavior 
would be for the damage coefficient to decrease as one- goes to very 
low energies. This could explain the low damage rate measured by 
Stabler and Curtin (ref. 23) for Q.27*MaV protons. Laboratory 
measurement of low energy proton damage, coupled with their analysis 
via the computer program,, is recommended. 

3. Extensions of the model to cover present developments appear 
feasible and should be undertaken. For example, a concentration gra- 
dient of lithium in lithium - doped solar cells (ref. 24) and the 
redegradation process appear to result in a non-uniform residual 
damage, which subroutine BOOT is uniquely able to handle for these • 
cells. In a separate extension, drift field solar cell performance 
In an AMO spectrum could be studied. 

4. The model should be applied to the evaluation of solar cell 
designs for specific flight raissiona, to relate design parameters 

to required performance -in space. Existing computer codes being 
used to fill this need are known to have deficiencies: e.g. the use 
of early and incorrect proton damage coefficients at intermediate 
energies, the inaccurate replacement of non-uniform damage effect by 
an "equivalent 1-M.eV electron" damage which is equivalent only -for 
a specific cell and a specific spectrum, the tedious calculation of 
covers lide shielding by Monte Carlo techniques, the neglect of 
changes in the junction parameters, etc. Such deficiencies in these 
earlier codes make use of the present <?ne advisable, although its 
results are no more precis® '$b§3J the current atata-of-ar t « 

' t » 

C. Kew Technology 

After a diligent review of the work performed under this 
contract, it was determined that no new innovation, discovery, 
improvement or invention. was develop*®* 
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GLOSSARY 


English Letters 

velocity of light in a vacuum (cm/sec) 

atomic diffusion coefficient of 'dopant atoms into silicon 
(cm 3 / sec) 

diffusion coefficient of electrons in p~type silicon . 

(cm 8 / sac)' 

diffusion coefficient of holes in n-tygs silicon (cm 2 /sec) 

electric field in cell due to impurity gradient (volts/ cm) 
energy of a particle in a solar cell 5 (MeV) 

.energy gap of silicon (1.11 sV) 

maximum residual energy of a proton after penetrating -a 
-thickness .t. (MeV) 5 • 

energy of a particle incident on a 3olar cell assembly 
(MeV) ..... , ■ 

rate of production of minority carrier per cm at depth x 
in silicon (carriers/cm 3 ~ sec) 

Planck r s constant (joule-sec or eV-sec) 

increment thickness (cm) 

spectral ir radiance outside solar cell assembly, at wave- 
length X (watte/ cm 2 - micron) 

Solar cell current through an external load '(milliamperes) 

Diode current -(mil Hamper es) ' 

photovoltaic current induced in the' solar cell (mllliamperes) 
current from a short circuited solar cell (milliamperes) 
diode reverse saturation current of a solar cell (milliamperes) 

minority carrier currant density (milliamperes /cm 2 ) 

damage coefficient relating proton or electron fluence 
to - degradation in L (dimensionless) 

Boltzmann 1 # constant (1.3B x 10 3 joule/molecule - K°) 

*f|a&d region minority carrier diffusion length (cm) 



L bass region minority carrier diffusion length before 
irradiation (cm) 

N concentration of dopant atoms (cm “) 
n parameter from Table 2 for proton range 

* *»3 

n electron concentration in p~type silicon (cm ) 
p hole concentration in n-type silicon (cm" 3 ) 

}?max maximum power available from a solar cell (milliwatts) 

q unit electric charge (1«602 x 10 19 coulomb) 

K range of a proton in silicon (mg/ cm®) 

R solar cell 3eries resistance (ohms) 

parameter from Table 2 for proton range 

T absolute temperature (degrees Kelvin) 

t process time of dopant diffusion to form n type surface 
region (sec) 

t covers lide thickness (cm) 

U light intensity incident on the solar cell assembly 
(milliwatts/cm 3 ) 

V solar cell potential across an external load (millivolts) 

V potential across an open circuited solar cell (millivolts)' 
oc 

V solar cell characteristic voltage (millivolts) 

x depth into the solar cell (cm) 

Xj depth of junction below solar cell surf ace. (cm) 

z depth into the solar cell assembly (cm) (z » t + x) 
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Greek Letters 


a(k) Absorption coefficient in silicon for light of wavelength 
X (cm*" 1 ) 

6 angle between a perpendicular to the solar cell surface 

and the direction of motion of the particle or photon 
being considered (degrees) 

9' angle of light ray in silicon, having angle 9 in space 

with respect to normal (degrees) 

X wavelength of light (microns) 

jx minority carrier mobility (cm 3 /volt - sec) 

T mean lifetime- of a minority carrier in the conduction 

band (sec) 

\ 

proton or electron fluonce to which the solar cell 
assembly has been exposed (cm* 2 ) 

Q resistivity of solar cell base region (ohm-cm) 

w 2tt times frequency of incident photon (sec”* 1 ) 


2 _ 'a 

P e 
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APPENDIX: . FORTRAN PROGRAM LISTING 

The computer program for solar cell performance predictions is 
listed here. In Fortran IV, it conforms to USA Standards .Institute 
standards. With, the' program entered in the computer,’ input variables 
corresponding to a particular solar cell' under consideration are to 
be entered by changing lines -150 to 270, 290, 940 to 970, and 990. 

The variables are identified in Table 3, on pp. 28 of this report. • 
Values greater than 200 for M, or greater than 50 for INA are not 
allowed in the program. 

A sample output is presented and discussed in Table 4, on pp. 33 
of this report. 
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lpo DIMENSION F.PC2Q0>*XL(2Q0)*HXC20Q)*DXC200)*C(20Q) 
110 DIMENSION 0 C200.) * FS C 50 ) * PI C 50 ) > OAC 50 > 

120 DIMENSION FOPC 1 0 )* PHP CIO)* EGEC 10)* PHEC 10) 

130 common/absc/dx/interv/hx/akt/c/gabsc/f 

140 COKM0N/DAMP/EP/DIFF/XL/ENEP./ES/FL1J/PI/AV1S/CA 

150 DATA XL3* D*RH0/1 50.*35.*10«/ ' 

160 DATA VO *p/43.*.l/ 

170 DATA T* XJU*CT/1 4» * .5*6./ 

1 SO DATA TEMP* U* THETA/300 . * I'A O'. jO.O/ 

190 DATA CE0PC I 5*1=1* 51/0.0*0.0*0.0*0.0*0.0/ 

200 DATA C ESP C I)* 1 = 6* 10 ) /O .' 0* 0 . 0* 0 . 0 * 0 .0* 0 . 0/ 

210 DATACPHPCI5* 1=1* 5) /O .0* 0 . Q* 0. 0* 0 . 0* 0 . 6/ 

220 DATA C PHP ( I ) * 1 = 6* 10 > /C . 0* 0 . 0* 0 . 6* 0 .0* 0 . 0/ 

230 DATACF9ECI )* 1=1* 5>/0.0* 0.0* 0.0*0.0* 0.0/ 

240 DATACECF.(I>*I = 6* 10 ) /O. 0 , Q . 0* 0 . 0* 0 . 0-* 0 . 0/ 

2 50 DATACPHECI)* 1=1* 55/0.0*0.0*0.0*0.0*0.0/^ 

260 DATACPHECI)* 1=6* ! 0) /0 . 0.- 0 . 0* 0 . 0* 0 .0* 0 . 07 
270 M-200 

2S0 n=M-2 
290 I MA= 50 
300 C 0 V= 1 • 0 

310 D=D*CTEMP/3Q0. )**<-! .5) 

320 XL0=XLS* 1 E- A’V C ?E*'P/3Q0 • 5**2. 

330 XJ!J=XJU*1 ,S~4 . 

340 CT=CT*2.54E~3 

350 T=T*2,54E-3 

360 IFCCT .NE« 0 * 0>'G0 TO 5 

370 C0V=0.0 

3S0 CT=XJU 

390 5 DC 10 I = l*M 

400 10 XL<n=XLC 

410 CALL ABSCIS CM,T*XJU5 

420 CALL LIGHT CM* TEMP* U* THETA) 

430 D0 30 0=1*10 
440 E3=E0PCIN5 
450 PHI = PHP C IN) 

460 I F ( E0 • L£ • 9 £~ 4 5 CO TS 30 

470 I F C EO *6E .200 •' 5 C3' TC 30 

480 IFCPHI .EO.O.O5G0 T0 30 

490 CALL COVER C E0* CT* PHI * XJU* C0V*-I NA> 

500 D0 20 MM=1*INA 

550 IFCESCMM5.LE..9E-45G0 T0 30 

520 CALL PROTON CXJU*M*MM* COV) 

530 CALL DAMAGE CM*MM*RH0) 

540 20 CONTINUE 
550 30 CONTINUE 
560 D0 60 IN= l * 10 
570 E0=E0ECIN) 

580 PHI=PHECIN5 

590 IFCE3.LE. .25)60 T0 60 

600 IF CEO .6E.4O.5G0 T0 60 

610 IFCPHI .EQ.-0.05 60 T0 60 

620 CALL ELECT < E0* PHI *M* CSV* XJU* CT* RH0 5 

630 60 CONTINUE ' 

640 C0-N=XLC2) 

6 50 CALL ROOT C N* M* D* XJU) 

660 CALL CURVE CD* TEMP* U* VO* R, CON* XL0* XJU* THETA) 

670 -STOP 
6S0 END' 
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690 * 

700 * 

710 * 

720 SUBROUTINE ABSGIS CM>T*XJ!J) 

730 DIKENSI0N DX< 200 ) *RX C200 ) 

7/‘0 C3KM0N /ABSC/DX/INTERV/HX 
750 FL0=FL3AT C M-20) 

760 DELTA=CT-XJU>/CFL0*CFL0-1 .)/2.*20. 5 

770 KK=M-1 

780 DX< I )=XJU 

790 -DO 10 I = I>20 

800 HX< I )=DELTA 

810 10 DXCI+1 )=DXCI)+HXCI) 

820 DO 5 1=2 1 * KK 

830 HXCn=FL0ATCI-2O>*DELTA 

sao 5 DXci+n=DXcn+Hxcn 

850 DX<M)=T 

860 RETURN 

870 END 

880 * 

890 * 

900 -* 
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910 SUBROUTINE: LIGHT CM* TEMP* U* THETA) . . 

920 DIMENSION A Cl 5 > * HC 1 5) * ST 1 5) * R <1 5 >* DK C 2001 * GC200 ) 

930 C3MM0N/ABSG/DX/GASSO/C 

940 DATA C A (I)*I=1*1 51/75000.V25800.* 11 GOO .* 7000 .* 4650* * 3330 . *2420 .* 
950+ 1690.* 1 120. *79 5. *380.* 1 80 . * 73 . * 20 .8* 4 . 4/ 

960 DATA CHC I >* 1 = 1 * 1 5)/ . 1 54*’ . 22* . 1 98* , 1 95* . 18 1 * * 1.62*. 144* 

970+ .127* . 1 12 7* .1003* .089 5* .0803* .0725* .0 665* .0606/ ■ 

980 IFCTHETA.EQ.O.O) G2 TO 20 
990 C I N= 3 . 8 

1000 AEG=THETA*.017453 
1010 AC=CSS (ARG5 
1020 AS=SINCARP) 

1030 AP=ATAN<AS/CCIN*<1 . - C AS/C IN ) **2 . 5 > ) 

1040 AS=C0S CAP) • 

10 50 G3 T0 30 
1060 20 AC=1 . 

1070 AS= 1 . 

1080 30 D3 10 1=1*15 
1090 H C I )=HCI )*U*AC/1 40 * 

1100 10 A<I>=ACI>*C3.323*EXPC-702.-6/TEMP> + l . >/l *319 
1110 P=6. 6258-34 • ' 

1120 V=2 . 998E8 
1130 DEN=P*V 
1140 W-0.4 
1150 1 = 1 

1160 C0E-O .05/3 . 

1170 SCI)=C0E*W*HCI>*A<I) 

1180 C0F.=C0E*2. 

1 190 3 5 l. ! = W+0 .05 
1200 1 = 1+1 

1210 S<I)=C0E*< 1 .+FL0AT<M3D(I+1*2> ) )*U*H<I7*ACX > 

1220 IFCI.LT* 14)60 T0 35 
1230 We W+0. 05 

1240 SO5!> = CC0E/2.>*W*H< 15>*A< 15) 

125OD0 45 J=!,M • 

1260 GG=0.0 
1270 D0 54 1=1*15 

1280 .R.CI).= $f I>*1 « E- 6#EXP C -AC I >*-0X( J) /AS) 

1290 54 GG=GG+R C I ) 1 

1300 45 GCJ5=GG/DEN 
1310 RETURN 
1.320 END 
1330 * 

13 40 * 

1350 * 
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1360 SUBROUTINE COVER CE0*CT*PHI*XJU*C0V*INA)' 

1370 DIMENSION ECS) *Ri C4)i ET1C4)* ESC50)*PIC50)*;GAC50) 
1380 DIMENSION ' ETC 4) * ROC 4> *R 1 C C 4) * ROC C 4) 

1390 C0MM0N/ENER/ES/FLU/P I /AN 6/CA 

1400 QATACE(I>*I=1*5>/0*Q* .3* *8* g.,*200./- 

1410 DATA C ROC ( I ) * I = 1 * 4) /2 . 8 1 * 3 . 945* 4 . 1 1 * 3 . 42/' 

1 420 DATA (ETC I )* 1= 1 *\4>/«99S* 1 .277* 1 . 46* 1 • 726/ 

1430 DATA CRI C< 15*1=1*4) /2 «07* 3*04* 3 « 18*2*85/. 

1 440 DATACETl (I)*I = t*4)/l«057> 1.375* ! .;576*1 .73/ ’ 

1450 IFCC0V* NE.O » 0) G0 TO 5 
1460 DO 6 1=1*4' 

1470 R1CI)=R0CCI)/2.33E3 
1480 6 ET1<I)=ETCI> 

1490 GO TO 7 ' 

1500 5 DO 47 1=1*4 

1510 47 R1CI)=R1CCI)/2.2E3 

1520 7 DO 3 1=1*4 . 

1530 IFCECI+1 5.GT.E05G0 TO 4 
1540 3 CONTINUE 
1550 GO TO 16 

1560 4 R0=R1 ( I)*CE0**£T1 CD) 

1570 fFCR0.LE.CT) GO T0 16 

1 580 £C=CE0**ETI CX')-CT/Rl CI))**U ;/ETlCI)> 

1590 IFCEC.LE..9E-4)60 TO 16 
1600 FAC=E0**E?1 Cl) 

36 10 AR=CT/CRlCn*CFAC-.9E-4**ETl CIO )> 

1620 TAN* SORT C I .-AR**2)/AR 
1630 AN=ATANCTAN> 

1640 STEP=AN/FL0ATCIWA> 

1650 XI 1 = 1 »/CFAC“EC*^ETl CD) 

1660 AR=0«0 

1670 D0 18 K=1*INA 

1680 CACK)=C0SCAR+.5*STEP) 

1690 ESCK>=CFAC-Ct/CR1 CI)*CACK) > )**U i/ETl < I ) > 

1700 XI2=XI1 

1710 XIl=FAC-C7/CRi CI>*C0SCAR-i-'STEP)) 

1720.XI1 = I ./CFAC-XIO 

1730 PICK>=PHl*CT*CXI-2-XI 1 )/Rl Cl) 

1740 18 AR=AR+STEP 
1750 XI1=*9£-4**ET1CI) 

1760 ^Xll^l ./‘CFAC-XI 1 .) 

1770 PI CINA5=PHI*CT J H(XI2~XI 1 >/Rl C I ) 

1780 G0 T0 19 
1790 ’16 ESC1)=0«G 
1800 19 RETURN ' 

1810 END 


1820 * 
1830 * 
1840. * 


Coining 7940: 
dsaaf £? of 3^98 


1510 47 £1(I)=R1C(I)/3.98K3 
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1 8 50 SUBROUTINE PR0T0N <XJU*M* MM* C0V5 

1860 DIMENSION HXC200) * EP (200)* ESC 50)* CAC 50) 

1870 DIMENSION EC 5-3 * RO < 4>* ETC 4) * ROC < 4) 

1 880 C0MM0N/INTERV/HX/DAMP/EP/ENER/ES/ANC/CA 
1890 DATA CECI 5 *1=1 *55/0.0* .3* .8* 2-. *200./ 

1900 DATACR0CCX5* 1 = 1 *45/2.8 1*3.945* 4. 11*3* 48/. 

1910 DATACETCI)* I=l*4>/«995* P.277* 1.46* 1.726/ 

1920 E°(U = ES(MM) 

1920 D0 47 1=1*4 

1940 47 R0CI5=R0C<I5/2.33E3 

1950 X=1 

1960 DB 21 1=1*4 

1970 IFCECI-U 5 «GT«EP< 1 > ) C0 2.2 
1980 21 CONTINUE 

1990 22 R0=CACMM)*RO( X >*< EPC 1 5**ETC I 5 5 • 

2000 IFCC0V.E9. 0.0560 T0 50 
2010 IFCR0.LT. XJU5G0 TO 16 

2020 EPC1 D=(EPC 1)**F.TCI>-XJIJ/CR0CI >*CACMM>) 5** Cl ./ETC I ) ) 

2030 50 IFCEPC15 .LE . 9E- 4) 60 T0 16 

2040 D0 6 K=2*M - 

2050 R0=CA<MM)*ROCI)*CEPCK-15**ET.CI> ) 

2060 IFCHXCK- 1) . GT.RSJGO T0 1 6 

2070 EP(K) = CEPCK-1 )**ETCI ) -HXCK- 1 5 / CRO C I ) *CA(MM9 ) •>**( 1 . /ETC I ) > 
2080 IFCEPCK5.LE..9E-45G0 T0 16 
2090 IF(EP<K).GT.ECI))60 T0 6 
2100 1 = 1-1 
2110 6 CONTINUE 
2120 60 T3 14 
2130 16 EPCK)=0 .0 
2140 14 RETURN 
2 1 50 END 
2160 * 

2170 * 

2180 * 

2190 SUBROUTINE DAMAGE CM*MM*RH0> 

2200 DIMENSION EPC200 ) * XLC 200 >* PI C.50) 

2210 C3MM0N/D IFF/XL /DAMP/ EP/FL U/PI 
2220 0ME=RH0**C-.?55 
2230 D0 19 I = 1 * M 
2240 IFCEPCI) .LT.3.5G0 T0 20 t 
2250 19 EPCI5=1.2E-5*0ME*<EP<I5**C-.95) 

2260 G0 T0 14 
2270-20 QEA-EXPC- l »04> 

228O.-0ME=0ME»4-.92EfJ5- * 

229Q.- D0421, . • • v.- 

2300 -IFCE'P-CJ) .i&.l .4-60 T0 '22 ^ ; 

2310 21 EPCJ)=0ME*DEA*CCEPC J)/.962>**<-.85^ ) 

2320 60 T0 14 • ...‘S' * ■ 

2330 22 D0 23 I=J*M 
2340 IFCEPC I > «LE. *9£~4> GS~ T0 14 
'2350 23 ,EPC I > = 0ME*EXPC-1 • OS’&EP'CI 5 5 
2360* 14 D0 5 I = 1 * M 
2370 IFCEPCU .EG. 0.0560 TB 24 
2380 A=XLCD*XLCI) 

2390 5’ XLCI)=SGRTCA/CA*PXCMM)*EP<I5+1 .0)5 
2 400.24 RETURN 
2410 END' ‘ 
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2420 * 

2430 * . 

2440 * 

2450 SUBROUTINE ELECT CE0*PHJ *M,C0V.*Xjy.» CT^RKO) 
2460 DIMENSION EPC200>*HXC200>,XLC200) ' 

2470 C0MM0N/DAMP/EP/I NTER V/HX/D I FF/XL 

2480 0ME=CUO./RH0>**.5)*1E-1O 

2490 EXN=-1O«/CE0**1.» 55 

2500 IF<E^*L£.l «p>G0 T3 10 

2510 DC«0ME*<1 a2*AL0GC2«17*£0)>**2. 

2520 G0 T0 20 

2530 10 DC~0ME*C o67*AL0GC4.>*E0> )**2* 

2540 20 TH=C0V*CT*2.2*XJU*2.33 
2550 EP<n=DC*EXPCEXN*TH5 
2560 A=XL<!>**2* 

2570 XLC*n = 5QRTCA/<A*PHI*EP< 1 > + l • 5 > 

2580 D0 40 l-2*M 

2590 TH=TH+HX<I-n^2.33 

2600 EPCX)=DC*EXPCEXN*TH> 

2610 A=XL<U**2,. 

2620 40 XLC-I)=SORT(A/CA*PHI^EPCI>f-i „5) 

2630 RETURN 
2640 END 
2650 * ‘ 

2660 # 

2670. * 


&ote; Instruction 2540 applies ta.s&l&t &02s with forcing 7940 
covers lides. For - sapphire wit ft: & daiisi-t^'r'^2! 3/<i9S gcs/cst s this 
bee#m®£ 



2680 SUBROUTINE ROOT <N,M,D,XJU> 

2690 DIMENSION HX C200) , DXC200> , XL <200 ) 

2700 DIMENSION 0 C200> > G <200 > , EPC200 ) 

2710 COMMON/ ft E S C / 0 X / 1 N T E F. V / H X / D I F F / X L 
2720 COMMON/ GABSC /G /DAMP / EP /AMT /C 
2730 DO 119 K= 1 , M 
2 7 40 HSUM*HXCK)+HXCK+I> 

2750 QX<K>«1 ./CXLCK+l >**2. >+2./CHXCK+!)*HXCK> > 
2760 XL < K) = 2 * / < HX < K > *HS UK ) 

2770 EPCK>=2./CHXCK+1 >*HS!JM> 

2780 DXCK}=DXCK>/EPCK> 

2790 XL(K)=XL(K>/EPCK) 

2S00 EP<K)=CCK+1 >/CD*EPCfO ) 

2810 119 CONTINUE 
2820 C<1>=0.0 
2830 CCSisl.OES 
2340 C2H-0 • 0 
2850 C2L=0.0 
2860 K=0 
2870 5 CONTINUE 
2S?0 CALL TRAP CN*K> 

2890 IF<CCK+2>> 10,20,30 
2900 20 IFCM-K)! 0*40*10 
2910 10 C2L=C<2> 

2920 C C 2 )=C C-2 >=?1 0 » 

2930 IF<C2L*C2H5Sd, 5* 50 
2940 30 C2H»C<2> 

2950 C C 2 > = C < 2 > /l 0 . 

2960 IFCC2L*C2H>50*5,50 
2970 50 C<25 = < G2L+C2H) /2 . 

2980 TEST=<C2H-C2L)/C2H . 

2990 IFCTEST.LT. l.OE-6><r.0 TO 40 
3000 CALL TRAP CN*K> 

3030 IFCC<K+2)>60,7Q,8O 
3020 70 IF (M“K > 60, 40, 60 
3030 60 C2L=CC2> 

3040 GO T0 50 - ' - 

3050 SO C2H=CC2> 

2060 GO TO 50 
3070.40 RETURN 
3080 END 
3090 * 

3100 * 

3110 ** 

3120 SUBROUTINE TRAP <N,K> 

3130 DIMENSION C C2D0 ) , DX(2Q0) , EP C 200 ) , XL <200 > 
3140 C0MM0N/ABSC/DX/DIFF/XL/DAMP/EP/AMT/C ' ■ 
3150 D0 99 K= 1 , N 

3160 C<K+2>«CCK+n*DXCK)-CCK>*XLCK)-EPCK> 

3170 IFCCCK+2> .GT.l .E3OIG0 TO 50 
3180 IFCCCK+2)'.LT.O.O)G0 TO 50 
3190 99 Continue 

3200- 50 RETURN 
3210 END 
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3220 * 

3230 * 

3240 * 

3250 SUBROUTINE CURVE CD, TEHP, U, VO* R, CON, XLO,XJU, THETA ) 

3260 DIMENSION C C200 > , HXC200 >' 

3270 COMM0N/AMT/C/INTERV/HX 

32S0 10 F0RMATC3HIL=,F6.2,3H MA> 1 OX, 2HU=, E6 . 1 , 9H MV?/CM**2> 
3290 20 F0RMATC3HVO=,F6.2,3H MV* tGX, 5HTEMP=, FS . G, 7H DEG. K) 
3300 30 F0RMATC2HR=,F6.3, 4H OHM, 1 OX, 6HTHETA=, F6 . 2, 5H DEG.) 
33 2 0 50 FORMAT C 3HI 0=, El 0 • 3, 3K MA, //) 

3320 70 FORMAT CSX, 5HI CMA) , 4X, 5HVCMV) , 4X, 5HPCMI’.') , / ) 

3330 80 FORMAT C7X,F6.2,3X,F6»1, 3X, Fo. 2) 

3340 G1 = 1 . 502 E- 2 9 

3350 UN= U*COS C THET A* .01 7453) 

3360 CU=Q1*D*CC2>*1.E3/HXC1) 

3370 CS= ♦ 425'43*AL0GCX JU*I E45-H .0057- 1 .0079*EXPC-53530 .*C0N> 
3380 CS=CS*CTEMP+400. >*UN*I E-4 
3390 CU=CU+CS 

3400 C 1 0= C XL0/CON) * C 1 . 572*AL0G< UN > - 3 .'36 ) *1 E5 
3410 Cie=CI3*EXP<-6492. 75/TEMP) 

3^20 WRITEC9, 3 0)CU, U 
3430 WR I TEC 9, 20) VO, TEMP 
3440 WRITEC9,30)R, THETA 
3450 WRITEC9, 50) CIO 
3460 t-.’R I TEC 9, 70 ) 

3470 ST£P=CU/1 0 « 

3480 DO 40 I=-l,10 
3490 XI=STEP*FL0ATCI-1 ) 

3500 .V=VO*AL0GC CCU+C'I0-XI)/CI0)-XI*R 

3510 P=XI*V*lE-3 

3520 40 WRITEC9,80)XI, V,P 

3530 STEP=CCU“XI)/1 0. 

3540 DO 60 J= 1,10 
3550 XI=XI +STEP 

3560 V=V,O*AL0GC CCU+CI S-XI ) /CIS ) -XI*R 

3570 P=XI*V*l£-3 

3580 60 UR-I TE C 9 , 80 ) X I , V, P 

3590 RETURN 

3600 END 
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